###单细胞分析
setwd("E:\\Ovarian\\1.picture\\2.scRNA\\CONDUCTED")
library(Seurat)
library(SingleR)
library(dplyr)
library(tidyverse)
library(patchwork)
library(ggplot2)
library(harmony)
library(devtools)
library(tidydr)
library(ggunchull)
library(paletteer)
library(ggsci)
###anchor整合样本
x=list.files()
dir=paste0(x,"/")
names(dir) = x      
counts <- Read10X(data.dir =dir)
scRNA1 = CreateSeuratObject(counts,min.cells = 3, min.features = 200)
scRNA.data=counts
scRNA=scRNA1

#QC 质量控制
#计算血红蛋白hemoglobin基因比例
HB.ref <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
rowSums(scRNA.data[HB.ref,])
HB.genes <- intersect(HB.ref, rownames(scRNA))
HB.genes

#计算每一个特征集的百分比，血红蛋白，线粒体，核糖体
scRNA[["percent.HB"]] <- PercentageFeatureSet(scRNA, features = HB.genes)
scRNA[["percent.MT"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
scRNA[["percent.Ribosome"]] <- PercentageFeatureSet(scRNA, pattern = "^RP[SL]")

##
col.num <- length(levels(scRNA@active.ident))
#画图展示每一个特征集在每个细胞中的百分比
p1=VlnPlot(scRNA,
           features = c("nFeature_RNA"), 
           cols =ggsci::pal_igv("default")(51), 
           pt.size = 0, #不需要显示点，可以设置pt.size = 0
           ncol = 1)+theme(legend.position = "",
                           axis.text.x = element_blank(),
                           axis.ticks.x =element_blank() )
p2=VlnPlot(scRNA,
           features = c("nCount_RNA"), 
           cols =ggsci::pal_igv("default")(51), 
           pt.size = 0, #不需要显示点，可以设置pt.size = 0
           ncol = 1)+theme(legend.position = "",
                           axis.text.x = element_blank(),
                           axis.ticks.x =element_blank() )
# p3=VlnPlot(scRNA,
#            features = c("percent.HB"), 
#            cols =ggsci::pal_igv("default")(51), 
#            pt.size = 0, #不需要显示点，可以设置pt.size = 0
#            ncol = 1)+theme(legend.position = "",
#                            axis.text.x = element_blank(),
#                            axis.ticks.x =element_blank() )
p4=VlnPlot(scRNA,
           features = c("percent.MT"), 
           cols =ggsci::pal_igv("default")(51), 
           pt.size = 0, #不需要显示点，可以设置pt.size = 0
           ncol = 1)+theme(legend.position = "",
                           axis.text.x = element_blank(),
                           axis.ticks.x =element_blank() )
p5=VlnPlot(scRNA,
           features = c("percent.Ribosome"), 
           cols =ggsci::pal_igv("default")(51), 
           pt.size = 0, #不需要显示点，可以设置pt.size = 0
           ncol = 1)+theme(legend.position = "right",
                           axis.text.x = element_blank(),
                           axis.ticks.x =element_blank() )
pdf('1_feature_UMI_MT_HB_Ribosome.pdf', width = 12, height = 10)
p1+p2+p4+p5+plot_layout(ncol = 2)
###把图片画到画板上面
dev.off()
#画图展示特征之间的关联性
p1 <- FeatureScatter(scRNA, 
                     feature1 = "nCount_RNA", 
                     feature2 = "nFeature_RNA",
                     cols =ggsci::pal_igv("default")(51) )+theme(legend.position = "",
                                                                 axis.text.x = element_blank(),
                                                                 axis.ticks.x =element_blank() )
p2 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", 
                     feature2 = "percent.HB",
                     cols =ggsci::pal_igv("default")(51))+theme(legend.position = "",
                                                                axis.text.x = element_blank(),
                                                                axis.ticks.x =element_blank() )
p3 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", 
                     feature2 = "percent.MT",
                     cols =ggsci::pal_igv("default")(51))+theme(legend.position = "",
                                                                axis.text.x = element_blank(),
                                                                axis.ticks.x =element_blank() )
p4 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", 
                     feature2 = "percent.Ribosome",
                     cols =ggsci::pal_igv("default")(51))+theme(legend.position = "right",
                                                                axis.text.x = element_blank(),
                                                                axis.ticks.x =element_blank() )

pdf('1_feature_relationships.pdf', width = 8, height = 6)
p1 + p2 + p3 + p4+plot_layout(ncol = 2)
dev.off()

#过滤细胞
#保留检测基因数在200至7000之间的细胞
#保留线粒体基因表达量小于5%的细胞
scRNA <- subset(scRNA, 
                subset = nFeature_RNA > 200 & 
                  nFeature_RNA < 7000 &
                  percent.MT < 10&
                  percent.HB < 3&
                  nCount_RNA < 100000)


##75508
###分割seurat对象
##分割成list按照样本
s.list <- SplitObject(scRNA, split.by ="orig.ident")
scRNAlist=s.list

for (i in 1:length(scRNAlist)) {
  scRNAlist[[i]] <- NormalizeData(scRNAlist[[i]])
  scRNAlist[[i]] <- FindVariableFeatures(scRNAlist[[i]], selection.method = "vst",nfeatures = 3000)
}


scRNA.anchors <- FindIntegrationAnchors(object.list = scRNAlist,anchor.features = 2000)


scRNA1 <- IntegrateData(anchorset = scRNA.anchors)

save(scRNA1,file = "IntegrateData.Rdata")

DefaultAssay(scRNA1) <- "integrated"

scRNA1=ScaleData(scRNA1)

scRNA1 <- RunPCA(scRNA1, verbose = T)
# t-SNE and Clustering
ElbowPlot(scRNA1, ndims=50, reduction="pca") 

scRNA1 <- FindNeighbors(scRNA1, reduction = "pca", dims = 1:20)
scRNA1 <- FindClusters(scRNA1, resolution = 0.5)


scRNA1 <- RunUMAP(scRNA1, reduction = "pca", dims = 1:20)#zuodaozhele
scRNA1 <- RunTSNE(scRNA1,reduction = "pca", dims = 1:20)
colors=c('#313c63','#b42e20','#ebc03e','#377b4c',
         '#7bc7cd','#5d84a4','#4B4B5D',"#EC7232")
DimPlot(scRNA1, reduction = "tsne", group.by = "orig.ident",cols = colors)+
  scale_x_continuous(breaks=NULL)+##去除网格线
  scale_y_continuous(breaks=NULL)+##去除网格线
  theme_dr()+
  labs(title = "Patient")
ggsave("1.patient.tsne.pdf",height = 4,width =5.7)
meta=scRNA1@meta.data
#############################################################################
DefaultAssay(scRNA1) <- "RNA"
scRNA <- ScaleData(scRNA1)

#save(scRNA,file = "anchor.tu.scal.Rdata")
#查看细胞周期的影响
cc.genes
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA))
#细胞周期评分
g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match = rownames(scRNA))
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = rownames(scRNA))
scRNAa <- CellCycleScoring(object=scRNA,  g2m.features=g2m_genes,  s.features=s_genes)

#查看细胞周期基因对细胞聚类的影响
scRNAa <- RunPCA(scRNAa, features = c(s_genes, g2m_genes))
DimPlot(scRNAa, reduction = "pca", group.by = "Phase",cols = ggsci::pal_npg("nrc", alpha = 0.6)(9))+
  #xlab(NULL)+ylab(NULL)+
  scale_x_continuous(breaks=NULL)+##去除网格线
  scale_y_continuous(breaks=NULL)+##去除网格线
  theme_dr()##箭头主题
ggsave("2_cellcycle_pca.pdf", width = 4, height = 4)





#logfc.threshold = 0.25    两组之间倍数改变大于0.25
#min.pct=0.25   至少在该cluster 25%的细胞中表达该基因
scRNA.markers <- FindAllMarkers(scRNA, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.25, 
                                logfc.threshold = 0.25
)
save(scRNA.markers,file ="scRNA.markers.Rdata" )
#挑选每个细胞亚群中特意高表达的10个基因
top10 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC) 

#整理成表格，只显示基因名字
top10_table=unstack(top10, gene ~ cluster)
names(top10_table)=gsub("X","cluster",names(top10_table))
write.csv(file="7_top10_marker_genes.csv",top10_table,row.names=F)
##SingleR，跑一编作为参考
pbmc=scRNA
counts<-pbmc@assays$RNA@counts
clusters<-pbmc@meta.data$seurat_clusters
ann=pbmc@meta.data$orig.ident
#ref=get(load("ref_Human_all.RData"))
ref=celldex::HumanPrimaryCellAtlasData()
save(ref,file = "singleR.ref.Rdata")
singler=SingleR(test=counts, ref =ref,
                labels=ref$label.main, clusters = clusters)
clusterAnn=as.data.frame(singler)
clusterAnn=cbind(id=row.names(clusterAnn), clusterAnn)
clusterAnn=clusterAnn[,c("id", "labels")]
write.table(clusterAnn,file="clusterAnn.txt",quote=F,sep="\t", row.names=F)

##科研猫
##cluster_tsne
DimPlot(scRNA, 
        reduction = "tsne",
        group.by ="seurat_clusters",label = T,
        label.size = 5,label.color = "black")+
  scale_color_manual(values =pal_igv("default")(51)[-1])+
  scale_x_continuous(breaks=NULL)+##去除网格线
  scale_y_continuous(breaks=NULL)+##去除网格线
  theme_dr()+
  theme(legend.position = "none")+
  labs(title = "Clusters")

ggsave(filename = "cluster_tsne.pdf",height = 5,width = 5)

save(scRNA,file = "scRNA.rdata")
##手动注释添加细胞类型
##dotplot
###dotplot
library(scRNAtoolVis)
sc2=scRNA
gene=c("EPCAM","CDH1","KRT7","KRT19",
       "PECAM1","VWF",
       "PDGFRA","LUM","COL1A1","MYLK","ACTA2","PDGFRB",
       "LYZ","S100A12","CD1C","CD68",
       "KIT","GATA2","CPA3","MS4A2",
       "NKG7","KLRD1","KLRB1","NCR1",
       "MS4A1","CD79A","CD19","JCHAIN",
       "CD3E","CD8A","CCL5","IL7R",
       "FOXP3", "TNFRSF4", "IKZF2", "IL2RA"
       
)
cluster=c(rep("Epithelial",4),rep("Endothelial",2),
          rep("Fibroblasts",6),rep("Myeloid",4),
          rep("Mast cells",4),rep("NK cells",4),
          rep("B cells",4),
          rep("T cells",4),
          rep("Treg cells",4))
y=data.frame(gene,cluster)
y$cluster=factor(y$cluster,levels = unique(y$cluster))
jjDotPlot(object = sc2,
          markerGene = y,
          dot.col = c('white','#4A3549','#0F3051'), #xtree = T,
          rescale = T, 
          anno = T,#textRot=45,textSize=8,lwd=0.1,#x.text.angle = 45,
          plot.margin = c(6,1,1,1), #base_size = 8,  
          rescale.min = 0,
          rescale.max = 1)
ggsave("DotPlot.pdf",height = 12,width = 10)



celltype=read.table("clusterAnn.txt",sep="\t",header=T,check.names=F)

scRNA@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
  scRNA@meta.data[which(scRNA@meta.data$seurat_clusters == celltype$id[i]),'celltype'] <- celltype$celltype[i]}
unique(as.vector(scRNA@meta.data$celltype))
##挑选想要得细胞
Idents(scRNA)="celltype"
#save(scRNA,file = "scRNA.Rdata")

#celltype_tsne
colors=c('#313c63','#b42e20','#ebc03e','#377b4c',
         '#7bc7cd','#5d84a4','#4B4B5D',"#EC7232")
DimPlot(scRNA, group.by="celltype", label=T, label.size=5, 
        reduction='tsne',pt.size = 0.4,cols = colors)+
  scale_x_continuous(breaks=NULL)+##去除网格线
  scale_y_continuous(breaks=NULL)+##去除网格线
  theme_dr()+
  theme(legend.position = "none")+
  labs(title = "Cell type")
ggsave(filename = "5.celltype.pdf",height =5 ,width = 5)

save(scRNA,file = "scRNA.rdata")
###
scRNA.markers.Treg=scRNA.markers[scRNA.markers$cluster==9,]
scRNA.markers.Treg=scRNA.markers.Treg[scRNA.markers.Treg$avg_log2FC>1,]
scRNA.markers.Treg=scRNA.markers.Treg$gene
save(scRNA.markers.Treg,file = "scRNA.markers.Treg.rdata")##258gene


#AUCell富集分析，制作基因集合
##BiocManager::install("AUCell")
#setwd("E:\\scRNA_brca\\2.scRNA\\1.Quality_control")
library(AUCell)
library(ggplot2)
library(Seurat)
library(clusterProfiler)


sc2=scRNA

cells_rankings <- AUCell_buildRankings(sc2@assays$RNA@scale.data, plotStats=TRUE) 
geneSets=list(Treg=scRNA.markers.Treg)
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings,nCores =1, aucMaxRank=nrow(cells_rankings)*0.1)

#X=getAUC(cells_AUC)
#Y=getAUC(cells_AUC)[geneSet, ]

aucs <- as.numeric(getAUC(cells_AUC))
sc2$AUC <- aucs
df<- data.frame(sc2@meta.data, sc2@reductions$tsne@cell.embeddings)
colnames(df)
class_avg <- df %>%
  group_by(celltype) %>%
  summarise(
    tSNE_1 = median(tSNE_1),
    tSNE_2 = median(tSNE_2)
  )
p1=ggplot(df, aes(tSNE_1, tSNE_2))  +
  geom_point(aes(colour  = AUC),size=0.2) + viridis::scale_color_viridis(option="F") +
  ggrepel::geom_label_repel(aes(label =celltype),
                            data = class_avg,
                            size = 2,
                            label.size = 1,
                            segment.color = NA
  ) + theme_bw()+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),
                        legend.position = c(.09, .19),axis.text.x = element_blank(),axis.ticks.x = element_blank(),##去除刻度线
                        axis.text.y = element_blank(),axis.ticks.y = element_blank(),
                        plot.title = element_text(color = 'black', size = 20, hjust = 0.5),
                        legend.key.size = unit(10, "pt"))+
  labs(title = "Treg-AUCscore",x="",y="")
p1
ggsave("Treg-AUCscore.pdf",p1,width = 5,height = 5)



unique(quantile(sc2$AUC))
##分高低表达两组画图
met=sc2@meta.data
sc2$AUC_group=as.vector(ifelse(sc2$AUC>0.5,"high","low"))

##kym画auc分组图
##immungAUC_Group

p2=DimPlot(sc2, 
           reduction = "tsne",
           group.by ="AUC_group")+
  scale_color_manual(values =pal_nejm("default", alpha = 0.6)(8))+
  theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),
        legend.position = c(.01, .1),axis.text.x = element_blank(),axis.ticks.x = element_blank(),##去除刻度线
        axis.text.y = element_blank(),axis.ticks.y = element_blank(),
        plot.title = element_text(color = 'black', size = 20, hjust = 0.5),
        legend.key.size = unit(10, "pt"))+
  labs(title = "Treg-AUCgroup",x="",y="")
p2+p1
ggsave("Treg-AUCgroup.pdf",p2,height = 5,width = 5)
ggsave("Treg-AUC.pdf",p2+p1,height = 5,width = 9)
###
scRNA3=sc2

###kym_gsva
# rm(list = ls())

library(Seurat)
library(ggplot2)
library(clusterProfiler)
library(GSVA)

gmtfile <- "h.all.v7.5.1.symbols.gmt"

hallmark <- read.gmt(gmtfile)

hallmark$term <- gsub('HALLMARK_','',hallmark$term)

hallmark.list <- hallmark %>% split(.$term) %>% lapply( "[[", 2)

expr <- as.matrix(scRNA@assays$RNA@counts)

dim(expr)

es.matrix = gsva(expr, 
                 hallmark.list, 
                 kcdf="Poisson",
                 method="ssgsea", 
                 abs.ranking=T ,
                 parallel.sz=3)


es.matrix =as.data.frame(t(es.matrix))
save(es.matrix,file = "gsva.rdata")
#es.matrix=es.matrix1
meta <-sc2@meta.data
#es.matrix=t(es.matrix)
meta=cbind(meta,es.matrix)
sc2@meta.data=meta

#FeaturePlot(scRNA,features = c("TGF_BETA_SIGNALING"),reduction = "tsne")

scRNA=sc2
group_list <- data.frame(sample = rownames(es.matrix), 
                         group = scRNA$AUC_group)
head(group_list)
library(limma)
design <- model.matrix(~ 0 + factor(group_list$group))
colnames(design) <- levels(factor(group_list$group))
rownames(design) <- rownames(es.matrix)
design


contrast.matrix <- makeContrasts(high-low, levels = design)

fit <- lmFit(t(es.matrix), design)

fit2 <- contrasts.fit(fit, contrast.matrix)

fit2 <- eBayes(fit2)

x <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")

head(x)

write.csv(x, "gsva.csv", quote = F)

df <- data.frame(ID = rownames(x), score = x$t)

#按照score的值分组
cutoff <- 10
df$group <- as.vector(ifelse(df$score>cutoff,"1",
                             ifelse(abs(df$score)<=cutoff,"2","3")))

#按照score排序
sortdf <- df[order(df$score),]
sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)
head(sortdf)

ggplot(sortdf, aes(ID, 
                   score, 
                   fill = group)) + 
  geom_bar(stat = 'identity') + 
  coord_flip() + 
  scale_fill_manual(values = c('dodgerblue4', 
                               'snow3', 
                               'palegreen3'), 
                    guide = FALSE) + 
  geom_hline(yintercept = c(-cutoff,cutoff), 
             color="white",
             linetype = 2, #画虚线
             size = 0.3) + 
  geom_text(data = subset(df, score < 0),
            aes(x=ID, 
                y= 0.1, 
                label= ID, 
                color = group),#bar跟坐标轴间留出间隙
            size = 2.5, #字的大小
            hjust = 0 ) +  #字的对齐方式
  geom_text(data = subset(df, score > 0),
            aes(x=ID, 
                y= -0.1, 
                label=ID, 
                color = group),
            size = 2.5, 
            hjust = 1) +  
  scale_colour_manual(values = c("black","snow3","black"), guide = FALSE) +
  xlab("") +ylab("t value of GSVA score, S-AUC High \n versus S-AUC Low")+
  theme_bw() + #去除背景色
  theme(panel.grid =element_blank()) + #去除网格线
  theme(panel.border = element_rect(size = 0.6)) + #边框粗细
  theme(axis.line.y = element_blank(), 
        axis.ticks.y = element_blank(), 
        axis.text.y = element_blank()) #去除y轴

ggsave("gsva.pdf",height = 10,width = 7)

scRNA=sc2
#计算差异marker
Idents(scRNA)="AUC_group"
marker=FindMarkers(scRNA,
                   ident.1 = "high",
                   ident.2 = "low")
save(marker,file = "findmarker.rdata")

write.table(rownames(marker),file = "diff.marker.txt",sep="\t",row.names=F,col.names=F,quote=F)
